function result = pictureSplit(filename)
    % close all;
    % clear;
    % filepath = "test2.png";
    % [filepath,name,ext] = fileparts(filepath);
    % filename = strcat(filepath, "/", name);
    I_cam = imread(filename);
    %figure;
    %imshow(I_cam);
    %load TEMPLATEDATA
    
    %I = I_cam(10+(1:460),90+(1:460));
    I = I_cam;
    picture = rgb2gray(I);
    
    addnum = 20;
    Isize = size(picture);
    temp = ones(Isize(1)+addnum,Isize(2)+addnum)*255;
    temp(1+addnum/2:Isize(1)+addnum/2,1+addnum/2:Isize(2)+addnum/2) = picture;
    picture = uint8(temp);
    % figure;
    % imshow(picture)

    % figure;%         U = edge(t,'canny');
    %         %[H, T, R] = hough(I,'Theta',-10:0.1:10);
    %         [H, T, R] = hough(U);
    %         P = houghpeaks(H,1);
    %         lines = houghlines(U,T,R,P);
    %         for m=1:length(lines)
    %             xy = [lines(m).point1;lines(m).point2];
    %             t = t.*~xor(t,imfill(t,xy(:,1)));
    %         end
    
    % subplot(2,2,1);
    % imshow(picture);
    se = strel('rectangle', [3,3]);
    Ibot = imbothat(picture, se);
    % subplot(2,2,2);
    % imshow(Ibot);
    I = Ibot;


    %H = fspecial('laplacian', 0);
    %hp1 = imfilter(I, H, 'replicate');
    %p1 = hp1;
    %S = I;
    %I = p1;
    makebw = @(I) im2bw(I.data,0.2);
    I = I < 30;
    % subplot(2,2,3);
    % imshow(I);
    
    I = (picture/2 + uint8(I)*127*1.2);
    
    % subplot(2,2,4);
    % imshow(I);
    % subplot(2,2,3)
    makebw = @(I) im2bw(I.data,median(double(I.data(:)))/1.2/255);
    I = ~blockproc(I,[92 92],makebw);
    %figure;
    %imshow(I);
    %% Remove Noise
    I = bwareaopen(I,20);
    %imshow(I);
    %% Clear the border
    I = imclearborder(I);
    %imshow(I);
    
    
    %% Find the largest box
    %hold on;
    R = regionprops(I,'Area','BoundingBox','PixelList');
    NR = numel(R);
    %hold off;
    
    maxArea = 0;
    for k = 1:NR
        A(k) = prod(R(k).BoundingBox(3:4));
        if R(k).Area > maxArea
            maxArea = R(k).Area;
            kmax = k;
        end
    end
    
    
    BBmax = R(kmax).BoundingBox;
    DIAG1 = sum(R(kmax).PixelList,2);
    DIAG2 = diff(R(kmax).PixelList,[],2);
    
    [m,dUL] = min(DIAG1);    [m,dDR] = max(DIAG1);
    [m,dDL] = min(DIAG2);    [m,dUR] = max(DIAG2);
    
    pts = R(kmax).PixelList([dUL dDL dDR dUR dUL],:);
    % h_pts = plot(pts(:,1),pts(:,2),'m','linewidth',3);%HUA TU XIANG KUANG
    
    XYLIMS = [BBmax(1) + [0 BBmax(3)] BBmax(2) + [0 BBmax(4)]];
    
    T = cp2tform(pts(1:4,:),500*[0 0; 1 0; 1 1; 0 1],'projective');
    IT = imtransform(double(I),T,'XData',[1 500], 'YData',[1 500]);
    
    I = IT;
    I = im2bw(I);
    %hold on;
    R = regionprops(I,'Area','BoundingBox','PixelList');
    
    %imshow(I)
    NR = numel(R);
    
    maxArea = 0;
    for k = 1:NR
        A(k) = prod(R(k).BoundingBox(3:4));
        if R(k).Area > maxArea
            maxArea = R(k).Area;
            kmax = k;
        end
    end
    
    
    BBmax = R(kmax).BoundingBox;
    DIAG1 = sum(R(kmax).PixelList,2);
    DIAG2 = diff(R(kmax).PixelList,[],2);
    
    [m,dUL] = min(DIAG1);    [m,dDR] = max(DIAG1);
    [m,dDL] = min(DIAG2);    [m,dUR] = max(DIAG2);
    
    pts = R(kmax).PixelList([dUL dDL dDR dUR dUL],:);
    %h_pts = plot(pts(:,1),pts(:,2),'m','linewidth',3);%HUA TU XIANG KUANG
    
    XYLIMS = [BBmax(1) + [0 BBmax(3)] BBmax(2) + [0 BBmax(4)]];
    
    %% Identify objects inside the box 
    A_tmin = 30; % Bounds for the digit pixel area
    A_tmax = 3600;
    digitbox_minarea = 20; % Bounds for the digit bounding box area
    digitbox_maxarea = 60^2;
    tsize = 24;
    pictures = zeros(2*tsize+1,2*tsize+1,9, 9);
    % figure;
    % imshow(I);
    
    cutsize = 5;
    
    kgood = zeros(1,NR);
    Pnew = zeros(NR,2);
    for k = 1:NR
        if R(k).Area < A_tmax && A(k) > digitbox_minarea && A(k) < digitbox_maxarea ...
                && R(k).BoundingBox(3) < 50 && R(k).BoundingBox(4) < 50 ...
                && R(k).BoundingBox(3) > 1 && R(k).BoundingBox(4) > 1
            
            Pnew(k,:) = [R(k).BoundingBox(1)+R(k).BoundingBox(3)/2 R(k).BoundingBox(2)+R(k).BoundingBox(4)/2];
            
            
            if inpolygon(Pnew(k,1),Pnew(k,2),pts(:,1),pts(:,2))
                %h_digitcircles(k) = plot(Pnew(k,1),Pnew(k,2),'ro','markersize',24);
                try
                    t = I(fix(Pnew(k,2)-tsize):fix(Pnew(k,2)+tsize),fix(Pnew(k,1)-tsize):fix(Pnew(k,1)+tsize));
                catch
                    continue
                end
        
                xpos = fix(Pnew(k,2)*9 / 500) + 1;
                ypos = fix(Pnew(k,1)*9 / 500) + 1;
                if (abs(0.5 - abs((xpos - 1 - Pnew(k,2)*9 / 500))) > 0.30 || abs(0.5 - abs(ypos - 1 - Pnew(k,1)*9 / 500)) > 0.3)
                    continue
                end
                % disp([xpos, ypos]);
                s = size(t);
                pictures(1+cutsize:s(1)-cutsize,1+cutsize:s(2)-cutsize,xpos, ypos) = t(1+cutsize:s(1)-cutsize,1+cutsize:s(2)-cutsize);
%                     figure;
%                     imshow(pictures(:,:,xpos,ypos))
            end
            
        end
    end
    % figure;
    rsize = 30;
    
    Rawpictures = zeros(2*rsize+1,2*rsize+1,9, 9);
    twoBitPics = zeros(28, 28, 9, 9);
    zeroPic = imread('zeroPic.png');
    
    for j = 1:9
        for k = 1:9
            % subplot(10,9,j*9+k);
            % imshow(pictures(:,:,j,k));
            t = I(max(1, fix((j-1)*500/9)+1-10):min(fix((j-1)*500/9+2*rsize-10)+1,500),max(1,fix((k-1)*500/9)+1-10):min(fix((k-1)*500/9+2*rsize)+1-10, 500));
    %         U = edge(t,'canny');
    %         %[H, T, R] = hough(I,'Theta',-10:0.1:10);
    %         [H, T, R] = hough(U);
    %         P = houghpeaks(H,1);
    %         lines = houghlines(U,T,R,P);
    %         for m=1:length(lines)
    %             xy = [lines(m).point1;lines(m).point2];
    %             t = t.*~xor(t,imfill(t,xy(:,1)));
    %         end
            %imshow(t);
            s = size(t);
            Rawpictures(1+cutsize:s(1)-cutsize,1+cutsize:s(2)-cutsize,j,k) = t(1+cutsize:s(1)-cutsize,1+cutsize:s(2)-cutsize);
            %imshow(Rawpictures(:,:,j,k));
            
            if sum(pictures(:,:,j,k))==0
                twoBitPics(:,:,j,k) = zeroPic;
            else 
                twoBitPics(:,:,j,k) = imresize(pictures(:,:,j,k), [28 28]);
                thresh = graythresh(twoBitPics(:,:,j,k)); %自动确定阈�??
                twoBitPics(:,:,j,k) = im2bw(twoBitPics(:,:,j,k),thresh); %对图像二值化
            end
            
        end
    end
    
    result = twoBitPics;
    
    